/* 多目标（高级）链路分析工具箱
 * 用途：
 * 输入一个发射机的位置和一组接收机（相同）的坐标和角度，以及其方向图，计算信号/干扰功率，对比接收机的灵敏度，并计算EMI余量
 * 参数
 * 1. freq                    分析频点
 * 2. Transmitter
 *    Transmitter.Power       发射机功率
 *    Transmitter.Gain        发射机方向图
 *                Gain.theta  一维数组theta
 *                Gain.phi    一维数组phi
 *                Gain.value  二维数组（dBi值）
 *    Transmitter.Position    坐标(x,y,z)
 *    Transmitter.M3          坐标变换M3
 *
 * 2. Receiver
 *    Receiver[i].Gain        接收机方向图
 *    Receiver[i].Position    接收机位置
 *    Receiver[i].M3          坐标变换
 *    ******************************************
 *    Receiver[i].Temp        接收机温度
 *    Receiver[i].LL          接收机线损耗   LL
 *    Receiver[i].NF          接收机噪声系数 NF
 *    Receiver[i].SNR         接收机信噪比   SNR
 *    Receiver[i].Band        接收机通信带宽 B
 *
 * 3. Option
 *    Option.TransformType    变换类型：0 - M4 变换矩阵；1 - 坐标 + 角度
 *    Option.Power            默认功率
 *    Option.Temp             默认温度
 *    Option.LL               默认线损耗    LL
 *    Option.NF               默认噪声系数  NF
 *    Option.SNR              默认信噪比    SNR
 *    Option.Band             通信带宽/速率 B
 *    Option.UF               频率单位（Hz/KHz/MHz/GHz/THz）
 *    Option.UP               功率单位（W/mW）
 *    Option.UL               距离单位（m/km）
 */

#function(freq, Transmitter, Receiver, Option)
#import toolbox;

// 计算传播信号功率函数
local t = exec('''app:/scripting/library/emc/em_prop.mx''');
local calc_prop = t.calc_prop;

// 默认选项
local OptionDefault = table();
OptionDefault.TransformType = 0;
OptionDefault.Power = 1;
OptionDefault.Temp  = 300;
OptionDefault.LL    = 3.0;
OptionDefault.NF    = 5;
OptionDefault.SNR   = 20;
OptionDefault.Band  = 1000;
OptionDefault.UF    = 1e9;
OptionDefault.UP    = 1.0;
OptionDefault.UL    = 1.0; // 空间距离单位

local upname = @(up)
{
    if(up == 1.0e+6) return "MW";
    if(up == 1.0e+3) return "KW";
    if(up == 1.0)    return "W";
    if(up == 1.0e-3) return "mW";
    if(up == 1.0e-6) return "μW";
    if(up == 1.0e-9) return "nW";
    return string.format("unknown power unit (%g)", up);
};

if(is_nil(Option))
{
    Option = clone(OptionDefault);
}else
{
    Option = table.union(Option, OptionDefault);
}

// table.view(Option);

local Temp = Option.Temp;
local Band = Option.Band;
local SNR  = Option.SNR;
local NF   = Option.NF;
local LL   = Option.LL;


///////////////////////////////////////////////////

// 计算所有发射机 - 接收机的接收强度
local GfuncMem = table();
local makeGfunc = function(GainData)
{
    local name = string(GainData);
    if(is_nil(GfuncMem[name]))
    {
        // 按dB计算的数据
        local Gfunc_dB = math.interp2(GainData.theta, GainData.phi, GainData.value);
        local Gfunc = @(theta, phi)
        {
            local x = Gfunc_dB(theta, phi);
            // println("Gain %g dB -> %g", x, math.pow(10, x/10));
            return math.pow(10, x/10);
        };
        GfuncMem[name] = Gfunc;
        return Gfunc;
    }else
    {
        return GfuncMem[name];
    }
};

// 波长使用国际单位制
local lambda = math.c0 / (freq * Option.UF);

local G1 = makeGfunc(Transmitter.Gain);
local G2 = [];
// local ReceivedPower = zeros(numel(Receiver));
local ReceiverPosition = zeros(numel(Receiver), 3);
for(local i=0; i<numel(Receiver); ++i)
{
    // 收集坐标
    ReceiverPosition[i, 0] = Receiver[i].Position[0];
    ReceiverPosition[i, 1] = Receiver[i].Position[1];
    ReceiverPosition[i, 2] = Receiver[i].Position[2];

    // 计算接收功率
    G2[i] = makeGfunc(Receiver[i].Gain);
    /*
    local p = calc_prop(lambda, G1, G2[i],
        Transmitter.Position * Option.UL, Receiver[i].Position * Option.UL,
        Transmitter.M3, Receiver[i].M3);
    ReceivedPower[i] = math.dB(p);
    */
}

///////////////////////////////////////////////////

// 构造智能数据控件
local idata = table();

// idata.ReceivedPower = ReceivedPower;
idata.ReceiverPosition = ReceiverPosition;
idata.Transmitter = clone(Transmitter);
idata.Receiver = Receiver;
idata.freq = freq;
idata.option = Option;

local tb = toolbox.new_smartdata(idata);
tb.title = _hT("LinkAnalysis");
tb.window = "toolbox.frame";


// 第一个操作：总体分析
// 绘制所有的接收机接收到的发射源的信号功率
local op = tb.add_operation(_hT("TotalAnalysis"));
op.autoexec = true;

local current_power = 1.0;

if(!is_nil(Transmitter.Power))
{
    current_power = Transmitter.Power;
}else
{
    current_power = Option.Power;
}

op.add_param(
    table {
        name  : "group_1a";
        label : _hT("TransmitterSettings");
        type  : "extent";
        value : "group"
    }
);

op.add_param(
    table {
        name  : "Power";
        type  : "string";
        label : _hT("Power");
        value : current_power;
        unit  : upname(Option.UP)
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "Frequency";
        type  : "string";
        label : _hT("Frequency");
        value : freq;
        unit  : util.ufname(Option.UF);
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "group_1b";
        type  : "extent";
        value : "group.end"
    }
);

op.add_param(
    table {
        name  : "group_2a";
        label : _hT("ReceiverSettings");
        type  : "extent";
        value : "group"
    }
);

op.add_param(
    table {
        name  : "Temp";
        type  : "string";
        label : _hT("Temperature");
        value : Option.Temp;
        unit  : "K";
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "NF";
        type  : "string";
        label : _hT("Noise_Factor");
        value : Option.NF;
        unit  : "dB";
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "SNR";
        type  : "string";
        label : _hT("Signal_Noice_Ratio");
        value : Option.SNR;
        unit  : "dB";
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "Band";
        type  : "string";
        label : _hT("Band");
        value : Option.Band;
        unit  : "Hz"; // bps
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "LineLoss";
        type  : "string";
        label : _hT("LineLoss");
        value : Option.LL;
        unit  : "dB";
        limit : nil;
        delay : 0
    }
);

// 警告级别（EMI Margin < 次阈值时标注黄色）
op.add_param(
    table {
        name  : "Critical_Level";
        type  : "string";
        label : _hT("Critical_Level");
        value : "10";
        unit  : "dB";
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "group_2b";
        type  : "extent";
        value : "group.end"
    }
);



op.add_param(
    table {
        name  : "group_3a";
        label : _hT("TransmitterPlatformSettings");
        type  : "extent";
        value : "group"
    }
);

// 发射机运动
op.add_param(
    table {
        name  : "Offset";
        type  : "string";
        label : _hT("Offset [X, Y, Z]");
        value : "[0,0,0]";
        unit  : util.ulname(Option.UL);
        limit : nil;
        delay : 0
    }
);

// 偏航角
op.add_param(
    table {
        name  : "Yaw";
        type  : "string";
        label : _hT("Yaw_Angle");
        unit  : "DEG";
        value : "0";
        limit : nil;
        delay : 0
    }
);

// 俯仰角
op.add_param(
    table {
        name  : "Pitch";
        type  : "string";
        label : _hT("Pitch_Angle");
        unit  : "DEG";
        value : "0";
        limit : nil;
        delay : 0
    }
);

// 滚转角
op.add_param(
    table {
        name  : "Roll";
        type  : "string";
        label : _hT("Roll_Angle");
        unit  : "DEG";
        value : "0";
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "group_3b";
        type  : "extent";
        value : "group.end"
    }
);

op.add_param(
    table {
        name  : "group_4a";
        label : _hT("PlotSettings");
        type  : "extent";
        value : "group"
    }
);

// 在场景图中显示增益
op.add_param(
    table {
        name  : "ShowGain";
        label : _hT("ShowGain");
        type  : "checkbox";
        value : 1;
        limit : nil;
        delay : 0
    }
);


local XRange = map(@(r){return r.Position[0];}, Receiver);
XRange.push(Transmitter.Position[0]);
local XWidth = max(XRange) - min(XRange);
// 绘制方向图时最大增益对应的距离大小（默认1）
op.add_param(
    table {
        name  : "GainRange";
        label : _hT("GainRange");
        type  : "string";
        value : XWidth / 15;
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "GainCutOff";
        label : _hT("GainCutOff");
        type  : "string";
        value : "-20";
        limit : nil;
        delay : 0
    }
);

op.add_param(
    table {
        name  : "group_4b";
        type  : "extent";
        value : "group.end"
    }
);

op.do_execute = function(idata, pdata)
{
    local Transmitter = clone(idata.Transmitter);
    local Receiver = idata.Receiver;
    local ReceiverPosition = idata.ReceiverPosition;
    local NumReceiver = ReceiverPosition.size(0);
    // local ReceivedPower = idata.ReceivedPower;
    local option = idata.option;

    local ReceivedPower = zeros(ReceiverPosition.size(0), 1);

    local g = fig.new_figure();

    local f = eval1(pdata.Frequency);
    local lambda = math.c0 / (f * option.UF);

    local Power = eval1(pdata.Power); // 发射功率
    local T    = eval1(pdata.Temp);
    local NF   = eval1(pdata.NF);
    local SNR  = eval1(pdata.SNR);
    local Band = eval1(pdata.Band);
    local LL   = eval1(pdata.LineLoss);

    local Yaw   = eval1(pdata.Yaw);
    local Pitch = eval1(pdata.Pitch);
    local Roll  = eval1(pdata.Roll);

    local Offset = double(eval1(pdata.Offset));

    Transmitter.Position += Offset;
    
    local m3 = math.transform_matrix("rotate", [0,0,1], Yaw)
             * math.transform_matrix("rotate", [0,1,0], Pitch)
             * math.transform_matrix("rotate", [1,0,0], Roll);
    m3 = m3[0:2, 0:2];

    // 计算接收功率
    for(local i=0; i<NumReceiver; ++i)
    {
        local p = calc_prop(lambda, G1, G2[i],
            Transmitter.Position * option.UL, Receiver[i].Position * option.UL,
            m3 * Transmitter.M3, Receiver[i].M3);

        ReceivedPower[i] = math.dB(p);
    }
    ReceivedPower += math.dB(Power);

    // 计算接收机灵敏度(dBmW)
    local Sensitivity = math.dB(math.boltzmann * T * Band / option.UP)+ NF + SNR + LL;

    local EMI_Margin = ReceivedPower - Sensitivity;

    // 对目标通信链路进行分级
    local EMI_Margin_Critical = eval1(pdata.Critical_Level);
    local index_sufficient   = find(EMI_Margin > EMI_Margin_Critical);
    local index_critical     = find((EMI_Margin <= EMI_Margin_Critical) .* (EMI_Margin > 0));
    local index_insufficient = find(EMI_Margin <= 0);


    fig.grid(3, 1);
    // 绘制坐标位置
    fig.plot(Transmitter.Position[:,0], Transmitter.Position[:,1], "linestyle", 0, "markerstyle", 1, "markercolor", "k", "markersize", 3);
    fig.info(_dT("Transmitter"));
    fig.plot(ReceiverPosition[:, 0], ReceiverPosition[:, 1], "linestyle", 0, "markerstyle", 2, "markercolor", "b", "markersize", 3);
    fig.info(_dT("Receivers (Total)"));

    // 绘制增益（只考虑水平面）
    if(pdata.ShowGain != 0)
    {
        local GainCutOff = eval1(pdata.GainCutOff); // 增益的截断阈值
        local phi_s = 0:360;
        local gain_s = zeros(1, numel(phi_s));
        local gain_range = eval1(pdata.GainRange);
        for(local i=0; i<numel(phi_s); ++i)
        {
            local phi = phi_s[i];
            local k = double([cosd(phi); sind(phi); 0]);
            local k1 = math.mldivide(m3, k);
            local theta1, phi1, r1 = math.cart2sph_deg(k1[0], k1[1], k1[2]);
            gain_s[i] = math.dB(G1(theta1, phi1));
            // println("Phi = ", phi, "| Gain = ", gain_s[i], "dBi");
        }

        gain_s[find(gain_s<GainCutOff)] = GainCutOff;
        local gain_m = (gain_s - GainCutOff) ./ (max(gain_s) - GainCutOff)* gain_range;

        local gvals = zeros(numel(phi_s), 3);
        gvals[:, 0] = (gain_m .* cosd(phi_s))[:];
        gvals[:, 1] = (gain_m .* sind(phi_s))[:];

        fig.plot(gvals[:, 0] + Transmitter.Position[0], gvals[:, 1] + Transmitter.Position[1], 'linestyle', 1, 'linecolor', 'b', 'markerstyle', 0);
        fig.info("Gain(dB)");
    }


    // 绘制临界区间的点
    fig.plot(ReceiverPosition[index_sufficient, 0], ReceiverPosition[index_sufficient, 1], "linestyle", 0, "markerstyle", 2, "markercolor", "green", "markersize", 3);
    fig.info(_dT("Receivers (Sufficient)"));
    fig.plot(ReceiverPosition[index_critical, 0], ReceiverPosition[index_critical, 1], "linestyle", 0, "markerstyle", 2, "markercolor", "purple", "markersize", 3);
    fig.info(_dT("Receivers (Critical)"));
    fig.plot(ReceiverPosition[index_insufficient, 0], ReceiverPosition[index_insufficient, 1], "linestyle", 0, "markerstyle", 2, "markercolor", "red", "markersize", 3);

    fig.xlabel(string.format("X (%s)", util.ulname(option.UL)));
    fig.ylabel(string.format("Y (%s)", util.ulname(option.UL)));
    fig.info(_dT("Receivers (Insufficient)"));

    fig.axis_equal(true);

    fig.title("Scene");
    fig.end_coord("Scene Plot");

    
    // 绘制所有的接收机的计算结果
    fig.plot(1:numel(Receiver), ReceivedPower, "linestyle", 1, "linecolor", 'k', 'markerstyle', 1, 'markersize', 3, 'markercolor', 'k');
    fig.info("Rx Power");
    fig.plot([1, numel(Receiver)], [Sensitivity, Sensitivity], "linestyle", 1, "linecolor", 'r', 'markerstyle', 0, 'linewidth', 2);
    fig.info("Rx Sensitivity");
    fig.plot([1, numel(Receiver)], [Sensitivity + EMI_Margin_Critical, Sensitivity + EMI_Margin_Critical], "linestyle", 1, "linecolor", 'purple', 'markerstyle', 0, 'linewidth', 2);
    fig.info("Rx Critical Level");
    fig.xlabel("Receiver No.");
    fig.ylabel(string.format("Power (dB%s)", upname(option.UP)));
    fig.title("Power & Sensitivity");
    fig.end_coord("Power & Sensitivity");

    // 绘制 EMI 余量
    fig.plot(1:numel(Receiver), EMI_Margin, "linestyle", 1, "linecolor", 'k', 'linewidth', 2);
    fig.info("EMI Margin");
    fig.xlabel("Receiver No.");
    fig.ylabel(string.format("EMI Margin (dB)"));
    fig.title("EMI Margin");
    fig.end_coord("EMI Margin");

    return g;
};

return tb;